#include "random.h"

#include <iostream>
#include <limits>
using namespace std;

#define A 10
#define k 64
#define N_ITERATIONS 20000
#define alpha 0.2
#define beta 0.4

class PRNG {
public:
  PRNG(unsigned int seed)
    : rng(*(new rng01_t(rng_t(seed), gm_uniform_double_01_t())))
  {}

  ~PRNG() {
    delete &rng;
  }

  double operator()() {
    return rng();
  }
private:
  rng01_t &rng;
};

double chi2_test_uniform(vector<unsigned int> &sample, unsigned int n_classes, unsigned int N) {
  //statistica di test
  double th = (1.0 / n_classes) * N;
  double v = 0;
  for(unsigned int i=0; i<n_classes; ++i) {
    double d = (double)(sample[i]) - th;
    v += d * d;
  }
  return v / th;
}

double chi2_test_discrete(vector<unsigned int> &sample, unsigned int n_classes, unsigned int N, vector<double> &th) {
  //statistica di test
  double v = 0;
  for(unsigned int i=0; i<n_classes; ++i) {
    double d = (double)(sample[i]) - (th[i] * N);
    v += (d * d) / (th[i] * N);
  }
  return v;
}

double serial(PRNG &rng, unsigned int order) {
  unsigned long n_tuples = 1;
  for(unsigned int i=0; i<order; ++i)
    n_tuples *= k;
  vector<unsigned int> z(n_tuples, 0);
  vector<unsigned int> r(order);
  unsigned int N = A * n_tuples;
  for(unsigned int i=0; i<N; ++i) {
    //genera order numeri
    for(unsigned int j=0; j<order; ++j) {
      r[j] = (unsigned int)floor((double)(rng() * k));
    }
    //incrementa il contatore della tupla
    unsigned long offset = 0;
    unsigned int u = 1;
    for(unsigned int j=0; j<order-1; ++j) {
      offset += r[j] * u;
      u *= k;
    }
    offset += r[order-1] * u;
    z[offset]++;
  }
  return chi2_test_uniform(z, n_tuples, N);
}

double chi(PRNG &rng) {
  return serial(rng, 1);
}

double gap(PRNG &rng, double a, double b) {
  unsigned int max_length_gap = /*k - 1*/ 20;
  //theoretical probabilities
  double p = b - a;
  double p_ = 1.0 - p;
  double sum = 0;
  vector<double> th(max_length_gap + 1);
  for(unsigned int i=0; i<max_length_gap; ++i, p*=p_) {
    th[i] = p;
    sum += th[i];
  }
  th[max_length_gap] = 1 - sum;
  /*
  cerr << "theoretical: ";
  for(unsigned int i=0; i<=max_length_gap; ++i)
    cerr << th[i] << " ";
  cerr << endl;
  */
  //compute the number of gaps to be considered
  double min_p = th[max_length_gap];
  for(unsigned int i=0; i<max_length_gap; ++i)
    min_p = (std::min)(min_p, th[i]);
  unsigned int n_gaps = (double)A / min_p;
  //cerr << "n. gaps: " << n_gaps << endl;
  //seek to a number in [alpha, beta]
  while(true) {
    double r = rng();
    if(r >= a && r <= b)
      break;
  }
  //random gaps
  vector<unsigned int> z(max_length_gap + 1, 0);
  for(unsigned int n=0; n<n_gaps; ++n) {
    //cerr << "gap #" << n+1 << " / " << n_gaps << endl;
    unsigned int length = 0;
    while(true) {
      double r = rng();
      if(r >= a && r <= b)
	break;
      length++;
    }
    z[(std::min)(length, max_length_gap)]++;
  }
  return chi2_test_discrete(z, max_length_gap + 1, n_gaps, th);
}

int main() {
  unsigned int seed = 5;
  PRNG r1(seed), r2(seed), r3(seed), rg(seed);
  for(unsigned int i=0; i<N_ITERATIONS; ++i) {
    cout << chi(r1) << ", " << flush
	 << serial(r2, 2) << ", " << flush
	 << serial(r3, 3) << ", " << flush
	 << gap(rg, alpha, beta)
	 << endl;
    cerr << (((float)i / N_ITERATIONS) * 100) << " %" << "\r";
  }
  cerr << endl;
}
